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A NEW PHASE SPACE DENSITY EOR QUANTUM 
EXPECTATIONS 

JOHANNES KELLER, CAROLINE LASSER, AND TOMOKI OHSAWA 


Abstract. We introduce a new density for the representation of quantum 
states on phase space. It is constructed as a weighted difference of two smooth 
probability densities using the Husimi function and first-order Hermite spec¬ 
trograms. In contrast to the Wigner function, it is accessible by sampling 
strategies for positive densities. In the semiclassical regime, the new density 
allows to approximate expectation values to second order with respect to the 
high frequency parameter and is thus more accurate than the uncorrected 
Husimi function. As an application, we combine the new phase space density 
with Egorov’s theorem for the numerical simulation of time-evolved quantum 
expectations by an ensemble of classical trajectories. We present supporting 
numerical experiments in different settings and dimensions. 


1. Introduction 

The wave functions describing the nuclear part of a molecular quantum system 
are square integrable functions on with specihc properties. They are smooth 
functions, but highly oscillatory and the dimension d is large. The frequencies of 
oscillations are typically related to a small semiclassical parameter e > 0, which 
can be thought of as the square root of the ratio of the electronic versus the average 
nuclear mass for the molecular system of interest. One expects 



ipix) {—ieVx)'ip{x) dx = 0(1) 


as e —>■ 0 for most nuclear wave functions ip G Often the semiclassical 

analysis of a molecular quantum system requires a phase space representation of 
the nuclear wave function, the most popular being the Wigner function 

}V^(z) := (27re)-" [ Pp{q + l)iP{q - De'^'^’/^dy, z = {q,p) G 
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The Wigner function is a square integrable real-valued function on the phase space 
T*'^d ^ many striking properties as for example 


x|i/>(x)pdx= / qWip{z)dz, 


ip{x){-ie'Vx)ip{x)dx = / pW^{z)dz. 


However, in most cases the Wigner function is not a probability density on phase 
space, since it may attain negative values. 

A guiding example is provided by the superposition of Gaussian wave packets. 
The semiclassically scaled Gaussian wave packet Qzq with phase space center zq = 
(<io,Po) € is defined as 

(1.1) gzo(x) := (7re)-‘^'^^exp(-^lx - qol^ + jpo ■ (x - ^qo)) , xeR’^. 

It satisfies 

/ d:lgz„(x)l^ dx = qo and / ^(x) (-ieV^Jg^Jx) dx = pg. 
jR'i jR-i 

Its Wigner function is a nonnegative Gaussian function centered at the point zg. 
However, the Wigner function of the superposition 

'4’ = 9zi + gz2y Zi,Z2 

has three regions of localization as seen in the left panel of Figure There are two 
regions around the points zi and 22 , respectively, where the Wigner function has 
nonnegative Gaussian shape, whereas in between around the midpoint of Zi and Z 2 
there is an oscillatory region with pronounced negative values. 
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Figure 1. Contour plots of the Wigner function (left), the Husimi 
function (middle), and the density (right) for a Gaussian super¬ 
position Ip = gzi +gz 2 ■ We chose the phase space centers zi = (0,1), 
Z 2 = (1, —§)) and the semiclassical parameter e = 0.14. Negative 
values are indicated by blue color. 
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One way of obtaining a nonnegative phase space representation of a wave function 
is to convolve its Wigner function with another Wigner function. One then calls 
the nonnegative function 

G 

a spectrogram of tjj- A widely used spectrogram is the Husimi function of ij), 

which is the spectrogram of ip with </> being the Gaussian wave packet go centered 
at the phase space origin. For the superposition example, the smoothing of the 
convolution removes the oscillations and widens the Gaussian profiles around the 
centers zi and Z 2 ', see the middle panel of FigureHowever, the smoothing also 
destroys important exact relations satisfied by the Wigner function; see also [121 
§7]. Let a : —)■ M be a smooth function with an appropriate decay property. 

Then, 

[ dx = [ a{q)W^{z)dz, 

jR'i 

while 

f a{x)\’4){x)\^ dx = f a{q)'H^{z)dz + 0{e) 

jR'i Jm.2d 

as £ —>■ 0, where the error term depends on second and higher order derivatives of 
the function a. Hence only the first moment of the position density x i-A |'0(a::)p is 
exactly recovered by the Husimi function. 

Our aim is now to systematically construct a new phase space density that 
maintains some of the key properties of the Wigner function, while being amenable 
to sampling strategies for positive densities. We propose a proper reweighting of 
the Husimi function and the spectrograms obtained from the first order Hermite 
functions 

Ve,ix) := {Tre)~'^^‘^^J^Xjexp{-^\x\‘^) , x G j = l,...,d. 

We define a real-valued function K by 

d 

:= (1 + ^ . 

By construction, is the difference of two nonnegative functions, a scalar multiple 
of the Husimi function * yVgg on the one side, and half the sum of the Hermite 
spectrograms > • ■ • > other side. For example, a single 

Gaussian wave packet 4 ’ = 5zo centered at the point zq G results in 

= (1 + i) (27>'e)”‘^exp(-A|^;_ 

- ^k-^o|^(27r£)-‘^exp(-A|2_2;o|2) ^ 

which is the difference of two well-localized positive densities. For the superposition 
of two Gaussian wave packets ip = gzi + gz 2 j we obtain a density p^ characterized 
by two islands of positive values around the centers zi and Z 2 -, which are surrounded 
by a sea of negative values; see the contour plot in the right panel of Figure and 
the explicit formula in §5.2[ 
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The new phase space function fj,^ allows for the exact representation of the 
moments of ip up to order three in the following sense. If a : —>■ K is a polynomial 
of degree less than or equal to three, then 


a{x)\tj;{x)\'^ dx = / a{q)fi^{z) dz, 
'4’{x)a{—ie\7)tf}{x)dx= / a{p)ii^{z) dz. 


For arbitrary smooth functions a : —)■ K and the associated Weyl quantized 

operator op(a), the quantum expectation value (i/>, op(a)i/>)j ^2 of the observable 
op (a) is approximated as 


(1.2) / ^j){x) op{a)tlj{x) dx = / a{z)fj,^{z)dz + 0{e^) 

JR<>' JR^t* 

for £ —> 0, where the error term depends on fourth and higher order derivatives of 
the function a; see Theorem |3 . 2| later on. Phase space approximations of quantum 
expectations and Wigner functions play a central role in the analysis of quantum 
systems, in particular in the semiclassical regime; see [121 §IV] or mi §7.1]. The 
idea of combining different spectrograms can also be found in the time-frequency 
literature; see e.g. [20] and the references given therein. However, the goal of [20] 
is cross-entropy optimization within a chosen set of spectrograms and not the ap¬ 
proximation of Wigner functions or quantum expectations. 

The second order accuracy with respect to e in the expectation value approxima¬ 
tion suggests using the new density in the context of molecular quantum dynamics. 
We consider the time-dependent Schrodinger equation 

i£dtilj{t) = + ' 0 ( 0 ) =' 00 , 

with a smooth potential function V : > K as provided by the time-dependent 

Born-Oppenheimer approximation. Let us denote by —)■ the flow of 

the corresponding classical equations of motion 

q = p, p = -\JV{q). 

Then, by Egorov’s theorem, we have 

(1.3) (0(t),op(a)0(f))^2 = / {ao^t){z)yV^„{z)dz + 0{£^) 

jR2d 


as £ —>■ 0, where the error depends on third and higher order derivatives of the 
functions o o and V as well as the L^-norm of the initial wave function ipQ. 
The Egorov approximation is computationally advantageous, in particular in high 
dimensions, since it allows to simulate the time-evolution of quantum expectations 
by an ensemble of classical trajectories. Over decades, it has been widely used 
in the physical chemistry literature under the name linearized semiclassical initial 
value representation (LSC-IVR) or Wigner phase space method. 

Our new phase space density comes into play here, since the combination of the 
approximations in and (|1.3[) provides 


( 0 (t),op(a) 0 (f ))^2 = / (ao$t)(z)AJ^o(z)dz-bO(e^ 
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which can be read as a new method for the computation of time-evolved quan¬ 
tum expectations by initial sampling from a difference of nonnegative phase space 
distributions; see Theorem |4.1| and the numerical experiments in 

1.1. Outline. Our investigation proceeds along the following lines. In ^ we briefly 
review phase space distributions as the Wigner function, spectrograms, and the 
Husimi function. ^ derives the new phase space density and proves our main 
result, that is, the second order approximation of expectation values by the phase 
space integration with respect to the new density function. Q applies this result 
to the quantum propagation of expectation values. Then, several explicit formulas 
for the density are derived in The numerical experiments in ^illustrate the 
applicability of the new approach for the dynamics of molecular quantum systems 
in dimensions d = 1, d = 2, and d = 32. Appendix |A| presents a sampling strategy 
for the density /x^ via the Gamma distribution used for the numerical experiments, 
while Appendix [B] provides further computational details. 


2. Phase space distributions 


In this section we review different possibilities for representing a square inte- 
grable function ^ G via real-valued distributions on the classical phase 

space Considering functions with frequencies of the order 1/e for a 

small parameter 0 < e <C 1, we work with the e-rescaled Fourier transform 

:= (27re)-‘^/2 f ^(g)e-iP''?Adq, p G 

We also use the Heisenberg-Weyl operator in £-scaling: 


Definition 2.1. The Heisenberg-Weyl operator associated with a phase space point 
z = {q,p) G is defined as 

:= e*P'(— -q), 'i/G 


Among its many striking properties, the following two will be important for us 
later on. We have 

and 

= exp (-^H(zi,Z 2 )) zi,Z 2 S 

where H x —>■ M denotes the standard symplectic form on i.e., 

(2.1) Vt{zi,Z 2 )'■= Jz 2 = qiP 2 -Piq 2 with J = 


0 Id 
-Id 0 


All phase space distributions considered here turn the action of the Heisenberg- 
Weyl operator on a wave function into a phase space translation by z, which is 
often referred to as a covariance property; see, e.g., (2.6) below. 


Remark 2.2. The Gaussian wave packet (1.11 with its phase space center at zq G 
is obtained by applying the Heisenberg-Weyl operator to the Gaussian 

goix) := (7r£)"‘^/^exp(-A|a;|2) 


centered at the origin, i.e., gzg = Tz„go- 
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2.1. Wigner functions. We start our discussion with the celebrated Wigner func¬ 
tion and recapitulate some basic relations. 

Definition 2.3. The Wigner function of a function if G is defined as 

W^{z) := (27r£)-'' [ fp{q + - ^)e^yP/^dy, z = {q,p) G 

JR‘^ 

Wigner functions are continuous square-integrable functions on phase space; 
however, they need not be integrable. The marginals are the position and mo¬ 
mentum density of the state, respectively. With a proper interpretation of the 
possibly not absolutely convergent integrals this means 

f = |i/’((7)P, [ W^{q,p)dq=\Te'if{p)\‘^, 

jR-i jR'i 

and in particular 

/ W4z)dz = Uf. 

jR2rf 

Wigner transformation preserves orthogonality in the sense that 

[ >V^(z)yV’</,(2:)dz = (27re)“'^ |('0,</>)|, ip, (f G L'^{R‘^), 

jRSrf 

and it turns the action of the Heisenberg-Weyl operator into a phase space trans¬ 
lation, i.e., 

Wt,v- = 

which is an example of the covariance property alluded above. 

Moreover, given a Schwartz function a : —>■ K, one can use Wigner functions 

to express expectation values of Weyl quantized linear operators 

(2.2) (op(a)i/;)(q) = (27re)"‘^ [ a(2±2i,p)V>(y)ed«-2^)'P/=dydp 

via the weighted phase space integral 

(2.3) (V', op(a)V') = [ a{z)Wji,{z)dz. 

We note that the oscillatory integral formula \2.2\ can be extended to more general 
classes of symbols a : —>■ K with controlled growth properties at infinity; see for 

example §4], §2] or [SI §2]. 

2.2. Spectrograms. Except for Gaussian states, Wigner functions attain negative 
values (see [5S]), and thus cannot be treated as probability densities. For example, 
any odd function if G L'^{R^) satisfies 

n;^{0) =-(2716)-'^ / |V^{|)|2dp<0. 

One way to obtain nonnegative phase space representations of a quantum state is 
to convolve its Wigner function with another Wigner function. 

Definition 2.4. Let if G L^(M‘^) and </> G 5(K'^). Then, is called a 

spectrogram of if. 
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In time-frequency analysis, spectrograms are typically introduced as the modulus 
squared of a short-time Fourier transform (see, e.g., the introduction in 0) so 
that the representation via the convolution of two Wigner transforms is derived 
subsequently. Spectrograms also form a sub-class of Cohen’s class of phase space 
distributions [6l §3.2.1.]. They satisfy 

(2.4) (W^*W^)(z) = (27re)-‘’*|(r,0_,i/;)^ z S 

where (j)-{x) := 4’{—x) for x € see also [HI Proposition 1.99]. Thus, spec¬ 
trograms are nonnegative and smooth by construction. The integrability follows 
from \2A\ by the square integrability of general Fourier-Wigner transforms z i—> 
{Tz(l>,ip) with (j),tp G L^(K‘’*), see Proposition 1.42 in [9]. Normalization is pre¬ 
served according to 

(2.5) / = 

Spectrograms also inherit the covariance property from the Wigner function, i.e., 

(2.6) = zGR^‘^. 

A particular spectrogram is obtained by convolving with the Wigner function of a 
Gaussian wave packet, which we will discuss next. 


2.3. Husimi functions. The most commonly used nonnegative phase space dis¬ 
tribution is the Husimi function; see e.g. [2 §4.1]. We consider the Wigner function 

(z) = (^£)-z G 

of the Gaussian wave packet go centered at the phase space origin and define: 
Definition 2.5. The Husimi function oi f) G is defined as the spectrogram 

H.^iz) := {yv^ *yVgg){z) = [ yv’^(r(;) (7re)“‘^e“l"’“’"l''/®dr(;. 


The Husimi function of ip is the spectrogram (2.4) with (p being the Gaussian 
wave packet go^ and since go is even, we have 


'H^(z) = (27re) \{T^go,ip)\" 


z G 


y)2d 


Also, since go is normalized, i.e., Ijgoll = 1) (2.5) gives 

[ 'H^{z)dz = 


That is, for ip G L^(R‘^) with jji/ijj = 1, the Husimi function is a smooth probability 
density on phase space. 

Remark 2.6. The Husimi function of ip G L^(M‘’*) is the modulus squared of the 
so-called Fourier-Bros-Iagolnitzer (FBI) transform, which associates with ip the 
mapping 

z^ (27r£)-^/"(T,go,iA). 

In contrast to the Wigner function and the spectrograms, the FBI transform is a 
linear, albeit complex-valued phase space representation; see pn Ghapter 3]. 
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Integrating a Schwartz function a : 
obtain 


against the Husimi function, we 




a{z)'HTp{z)dz = [ a{z){W- 4 , *Wgf^)[z)dz = [ {a yVgg){z)yV^{z)dz 

jRSti jR2d 

= (-0,op(a* WgJV’), 


where the last equation uses the duality relation (2.3) between Weyl quantized 
operators and the Wigner transform. Therefore, 


/R2d 


a{z)n^{z)dz = ('(/',opa,^(a)V’), 


where 


oPaw(a) := op(a * Wgo), ae5(K ) 


denotes the anti-Wick quantized operator of the function a; see for example [HI §2.7] 
and [U §11.4]. Weyl and anti-Wick quantization are e-close in the following sense: 


Lemma 2.7. Let a : —>■ K &e o Schwartz function and e > 0. Then, there 

are two families of Schwartz functions rf,r| : —>■ M that depend on fourth and 

higher order derivatives of a, so that 

oPaw(a) = 0P(a + l^a) + op(rf), 

0Paw(a - |Aa) = op(a) -P op(r|), 
where sup^^g I! op(^PIIl 2 _>l 2 < oo for both j = 1, 2. 

Proof. The lemma is essentially proven in [181 Proposition 2.4.3] or [HI Lemma 1], 
and hence we only sketch the proof for the second of the two equivalent identities. 
We write out the definition 


0Paw(a - I Aa) = op {Wg„ *{a- |Aa)) 
and Taylor expand a — |Aa around z in the integral 

Wg„ * (a - |Aa) = i7Te)-‘^ [ (a - |Aa)(C)e-l^-f ^^dC. 

Due to the symmetry of the Gaussian, all Taylor expansion terms with odd deriva¬ 
tives of (a — |Ao) vanish. The computation 


l^iV<i(2a)! 



|Aa) 


implies the second order approximation 

Wffo * (a - |Aa) = (a - |Aa) -P |A(a - |Aa) -k 0{e^) 

= a + 0{£^), 

where the 0{£^) term is of Schwartz class. Applying the Calderon-Vaillancourt 
Theorem (see, e.g., [9l §2.5]) concludes the proof. □ 


Remark 2.8. The result of Lemma 2.7 can formally be read in terms of the heat 
semigroup {exp(tA)}t>o as 


where a : 


a * Wgg = exp(|A)a = a -I- |Aa -k 0{e^), 
is a Schwartz function. 
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3. The new phase space density 


We learn from the preceding discussion of phase space distributions, in particular 
from Lemma |2.7[ that the Husimi function allows to approximate an expectation 
value according to 

(V',op(a)'(/') = [ (a - lAa){z)'Hjij{z)dz + 0{e^) 

as £ —> 0, where the error depends on the fourth and higher order derivatives of a 
and the L^-norm of ip. An integration by parts provides 

{ip,op{a)ip) = [ a(z) - |A'H^)(z)dz + 0(£^), 

and motivates us to define the following new phase space density. 

Definition 3.1. For ip G we define the phase space density p,,/, : 

:= — |AH^. 

We summarize the key property of the new density as follows: 


Theorem 3.2. Let a : —>■ K fee a Schwartz function. Then, there exists a 

constant C > 0 depending on fourth and higher order derivatives of a such that for 
all Ip G 


{ip,op{a)ip) - / a(z)/x^(z)d 0 <Ce^ 

where the density —>■ M was defined in Definition 

Proof. By Lemma |2.7[ we have 


with 


3.1 


{ip,op{a)ip) = [ {a-lAa){z)H^{z)dz + s‘^{ip,op{r'^)ip) 


|(ife,op(r^)V')| < C\\ipf, 

where the constant C > 0 depends on the fourth and higher order derivatives of a. 
The Husimi function is smooth and bounded since 

= 

/ a{z)fj,.^{z)dz. 


n^{z) = (27r£)-^ \{g^,iP)\^ < i27re)-^\\g^\n 
Hence, integration by parts implies 

[ {a - jAa){z)'H^{z)dz = [ a{z){'H.^ - ^A'H^){z)d. 

Therefore, 


{ip,op{a)ip) - / a{z)yL.ii,{z)dz 


/R2‘i 


= e^|(ife,op(r^)i/>)| < 


□ 


is a polynomial of degree less than or equal to 


which concludes the proof. 

Remark 3.3. Whenever a : 
three, the constant C > 0 of Theorem |3.2| vanishes so that the phase space integra¬ 
tion with respect to exactly reproduces the expectation value. Moreover, using 
higher order Hermite spectrograms, one can also construct densities which yield 
approximations of expectation values with higher order errors in £, or, equivalently, 
which are exact for polynomial symbols a : K of higher degree. We refer 
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to the thesis m §10.5] for the next order result and an outline on how to prove 
higher order approximations. 

Our next aim is to derive an alternative expression for the new density showing 
that it is a linear combination of spectrograms. 

3.1. The new density in terms of Hermite functions. The Laplacian of the 
Husimi function can be related to the Wigner function of the £-rescaled first order 
Hermite functions as follows: 

Proposition 3.4. Let £ > 0, j £ {1,..., d}, and 

ipe^ix) := (7r£)"‘'/^y^a;jexp(-i|a;|^) , a; G K'', 
be the first order Hermite functions. Then, for all if G 


j=i 


and consequently 


I AHv, = (1 + i)H^ ■ 

j=i 

Proof. Let z G We compute 

AWgfiz) = (7r£)-‘^Ae-l"l'/" = (7r£)-^ V • (-fze-'"''/") 

= (-f+ ^|^|2) e-NV.. 

Moreover, by direct computation or [13 Theorem 1], 

W^^^.(z) = -(^£)-‘^(l-f|z,p)e-N^/^ 

such that 

E (2d - d - f 


i=i 


= -(7r£)-^ (2d-f|zp) e-l^lV- + d.W,„(z) 


and 


AW,„(z) = -f(^£)-‘^ (2d - fjzHe-l^l^/^ = - ¥ Wao(^)- 


j=i 

To conclude the proof we note that AH^ = A(>V^ * Wg^) = W,/, * AWg^,. □ 

Remark 3.5. For any fj G the real-valued density is a weighted sum of 

spectrograms and therefore smooth and integrable. It satisfies the normalization 
condition 

[ /x^(z)dz = llV’f 

and the covariance property 

/Jx’s'i/? — MV'(* z G M 

Next, we add a further characterization of the new density that does not explicitly 
require a convolution. 
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3.2. The new density in terms of ladder operators. The Gaussian wave 
packet go centered at the origin and the first order Hermite functions ,..., ^Ped 
can be characterized by the raising and lowering operators 

A^ = ^{x-e^d:) and A=^{x + eWa;), 

respectively. On the one hand, we have 

spanjgo} ={'*/'€ | Ajtlj = 0 for all j = 1,..., d} , 

for the kernel of the lowering operator A = (Ai,..., A^). On the other hand, the 
components of the raising operator applied to the Gaussian wave packet go generate 
the first order Hermite functions in the sense that 

‘Pej = A^jgo, j = l,...,d. 


Proposition 3.6. For all ip £ j = 1, ... ,d, and z = {q,p) £ we have 

|(TzA]go, 

= ( 27 re) (^g^,(^Aj - , 

where z^ := q + ip a and g^ is the Gaussian wave packet defined in (1.11. 
Consequently, 

di>{z) = (27re)"‘^ I (1 + ^)|(r^go,'0)|^ “ Ii/-)| 

\ 


( 27 re) + i)\{ 9 ^^'>P)\^ - \ '^\{ 9 z,(Aj - 


j=i 


Proof. The relation ( |2.4[ ) implies 

= (27r£)-" \{T,A]go,iP)\^ 

since {ipej)-{x) = Pefi—x) = —ipefix) for X £ For the second representation of 
the Hermite spectrogram we compute 


sdj o Tz = ipjT^ + Tz o edj 

and deduce 

= fke ~ ° 

= H - iPj)) ° 

so that 

{TzA]go,ip) = {A]g^,ip} - -^{qj + ipj){g,„ip}. 

□ 


Proposition 3.6 will be used for explicit expressions of later on in ^ when ip 
is a superposition of Gaussian wave packets or a higher order Hermite function. 
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4. Quantum dynamics 

As an application of the new density we consider the approximation of expecta¬ 
tion values for the solution of the time-dependent semiclassical Schrodinger equation 

iedtijj(t) = Htpit), V’(O) = i/'o, 

where the Schrodinger operator H = op{h) is the Weyl quantization of a smooth 
function h : —)■ K of subquadratic growth, that is, all derivatives of the func¬ 

tion h of order two and higher are bounded. Then, H is essentially self-adjoint 
(see [25l Exercise IV. 12]) so that for all square integrable initial data ipo G 
there is a unique global solution 

t G K. 


The classical counterpart to the Schrodinger equation is the Hamiltonian ordinary 
differential equation 


i(t) = jyh{z{t)) with J 


0 

-Id 


Id 

0 


^ j^2c/x2(i 


The associated Hamiltonian flow ^ is globally defined and smooth 

for all times t G K, since h is smooth and subquadratic. 

In this setup we obtain the following quasiclassical approximation of time evolved 
quantum expectations using the new phase space density. 


Corollary 4.1. Suppose h : ^ M is a smooth function of subquadratic growth 

and H = op{h). Let ip G with = 1- Then, for all Schwartz functions 

a : —?► M, and t G ffi, there exists a constant C = C{a, h,t) > 0 such that 


3 op(a)e - [ {a o ^t){z)^j,^{z)dz 


< Ce^ 


with the density from Definition 3.1 where —>■ is the Hamiltonian 

flow associated with h. 


Proof. The crucial element of our argument is Egorov’s theorem [31 Theorem 1.2], 
which provides 

e’"*/" op(a)e-‘^‘/" = op(a o + 0(6^), 

where the error depends on third and higher order derivatives of a o and h, 
respectively. This means for the expectation value 

op(a)e-‘"‘/^V’) = op(a o + 0(e^). 

Now it remains to apply Theorem |3.2| to obtain 

op(a)e-‘^*/^i/>\ = / (a o 4>t)(z)/r^(z)dz -H O(e^). 

' ! 7B2d 

□ 


Remark 4.2. Replacing the new density by the Husimi function in Corollary |4.1| 
deteriorates the approximation in the sense that 
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as £ —>■ 0. It requires an additional system of coupled ODEs involving higher order 
derivatives of the Hamilton function h to retain second order accuracy with respect 
to £; see m Theorem 3]. 


Remark 4.3. Since the Hamiltonian h is preserved by the classical flow the 
constant C{h,h,t) = C{h,h) of Corollary 4.1 does not depend on time, so that 
the approximation error of the total energy expectation value is of size 0(£^) but 
time-independent. 


Remark 4.4. In the special case of a harmonic oscillator h{z) = z’^Az, with A G 
j^ 2 dx 2 d positive definite, generating a flow 4>t that is a linear orthogonal map on 
phase space, one can easily see that 


In other words, satisfies the classical Liouville equation. In general, the time 
evolution of is much more intricate, see the equations for the Husimi function 
derived in [21 Theorem 4.5 and §4.2]. 


5. Examples of Phase Space Densities 


In this section we explicitly compute the new density from Definition 3.1 
in three different cases, namely when ip is a Gaussian wave packet, a Gaussian 
superposition, or a multivariate Hermite function. 


5.1. Gaussian wave packets. The Gaussian wave packet 

9z{.x) = (7r£)"^/'‘exp(-i|x - qp + '-p ■ {x - \q)) , xGR^. 

centered at z = {q,p) S has the Wigner function 

(5.1) = (7r£)-'^exp(-i|w-z|2) , w e 

Its Husimi function is a Gaussian function, too, but broader, that is, 

TLg^{w) = (27r£)“‘^exp(—— zp) , w G 

The Hermite spectrograms of can be related to another Husimi function. Indeed, 
by the covariance property of the spectrograms, 

(Ws. * = (Who * - ^) = 

= ( 27 re)-‘^i|wj -Zjf exp(-i|w;-zp) 

for all w = {wi, ..., Wd) G with wi,... ,Wd G and all j = 1, ..., d. Summing 
all the Hermite spectrograms then yields 

dgAw) = ( 27 re)"'^ (l + f - ik - A^) exp(-^|ic - z^) , w G 

Eigurej^ plots the new density /dg^(w) together with the corresponding Wigner and 
Husimi function in terms of the distance \w — z\. Its polynomial prefactor puts the 
density Pg^ in between the Wigner and the Husimi Gaussian. 
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Figure 2. Plots of the Wigner function (black dashed), the 
Husimi function (blue dashed), and the new density (red solid) 
for Gaussian wave packet t/j = in dimension d = 1 with semi- 
classical parameter e = 0.14. We plot the three functions in terms 
of the distance \w — z\. 

5.2. Gaussian superposition. Next we compute the new density for a Gauss¬ 
ian superposition, that is, for 

V' = 5zi+ffz2, 21 , 2:2 G 

Writing the value of the Husimi function at a point 2 G as 

n^.iz) = {27re)~'^\{g^,g^,) + \{gz,gz2)\‘^ 

= (27re) (|(< 7 zi 52 i)| + 1(52 > 522 ) I ~i~iidzi, gz){gzT gz 2 ))) 

motivates us to derive an explicit formula for the inner product of two Gaussian 
wave packets. 

Lemma 5.1. For any 21,22 G we have 

( 521 , 522 ) = exp + ^ 0 ( 21 , 22 )), 

where H is the standard symplectic form on T*M.‘^ = as defined in 
Proof. By a direct calculation one can show 

( 50 , 52 ) = {go^T.go) = exp (--*1^), ^ 

For the general case we then have 

— (^21 5 o 7 ^225^0) — (5^07^—21^22^0) 

= exp (^ 11 ( 21 , 22 )) ( 50 ,Tz2-2i5o) = exp (^^^( 21 , 22 ) - 


( 2 . 1 ). 


□ 
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By Lemma |5.1| we have 

{ 9 z^,gz){gz,gz 2 ) = exp exp (^0(zi - Z2,z)) 

and obtain for the Husimi function 

n^{z) = (27re)-‘^ (exp + exp 


2s 


+2exp(- l" "^1' ) cos(if2(zi -Z 2 ,z))) • 


For the Hermite spectrograms we first use Lemma |3.6| to obtain 

2 


iW^*W^){z) = i27rey 


fc = l 


F] ( j^O 


We notice that, for j = 1,..., d and m = 1,2, 




where := Qmj + ipmj & C. This implies that 

[gz, = (^g,, A,)go 

--z^ 

^ " ^9z,gzm)- 

Consequently, since \z^ j — z^\ = \zm,j — Zj\ for each m and j, we have 


(27r£)"(W^*W^^ J(z) = 


^9z,gzi) + {gz,gz2) 


Now observe that 


i=i 

Therefore, 

d 

i=i 

- 1 / CXD 

“ (27r£)‘i 1 2s 


= 


+ iRe 


= E[feu 

i=i 

= (z- Zl) 

-Qj)- Kpi 

• (z - Z2) +] 




lexp(- 


\z-Zi\'^ + \z-Z 2 \ 

4s 


•) (z - Zi) ■ (z - Z2) COS (Afi(z^ - Z2, z)) 
- fi(z - zi,z - Z2) sin (^^(zi - Z2, z)) 


Combining the Husimi function and the Hermite spectrograms gives the density 
as plotted before in Figure in the Introduction. 
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5.3. Hermite functions. Let us consider the higher order Hermite functions 
:= A; = (fci,...,fcd) e N"', 

which result from the A-fold application of the raising operator to the Gaussian 
wave packet go centered at the origin. The FBI transform of a Hermite function is 
known as 


(5.2) 

for z e K 

(5.3) 


(2^e)-'^/2(5„^,) = (^(g-ip))'exp(-i|zp) 

see la § 2 ] or m Proposition 5]. Hence, the Husimi function satisfies 

n J, ^ 




1 

(27re)‘^A! 




exp(-i|z|- 


= Y[hk^{zj), 


where 

hn{w) ■■=ic e 

denotes the Husimi function of the univariate n-th Hermite function n € N. 
For the multivariate Hermite spectrograms we obtain the following: 

Lemma 5.2. For all k and j = 1,d and z G we have 


= {kjhk.-i{zj) - 2kjhk^{zj) + {kj + l)hk-+iizj)) ■ Afc„(z„). 


Proof. Since Ajtpk = ^y~kj(fk-ej. Proposition 3.6 implies 

(Wva. * W^e^)(z) = (27re)-"* {g^,Ajipk) - ^zf{g^,ipk) 

= {2tT£) \/^{gzi Pk — ej) ~ {9ztFk) 


By the formula (|5.2|) for the FBI transform, we obtain 

k — Er. 


(W^.*W_)(^) = 




(27re)^fc! 

Then it remains to observe that 


exp(-d-|z|2) 

(27re)^fc! 


- 

1 






^z- 




2ki-2 


k—e 


(27re)Aj! 

= kjhkj—i{zj^ 2kjhkj{zj) {kj T l')hk^-\-i{zj'). 


□ 


The com bination of the Husimi function (5.3) and the Hermite spectrograms of 
gives an explicit formula for the density when if = tpk for k G 


Lemma 


5.2 
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6. Numerical Experiments 


We present numerical experiment^ for computing expectation values for the 
solution of the time-dependent Schrddinger equation 

iedtipit) = + •0(0)= V'o, 

with three different potentials V : —)■ K. The various setups shall illustrate 

important aspects of our new algorithm, such as the second order accuracy with 
respect to s, the good applicability in higher dimensions, and the capability of 
describing fundamental quantum effects. 

6.1. Discretization. For the algorithmic discretization of Corollary |4.1| we proceed 
similarly as in na nni [H]- We consider various smooth functions a : —>■ M 

and evaluate the phase space integral on the right hand side of the semiclassical 
approximation 

(0(t),op(a)0(t)) = [ {ao^t){z)iJ.^g{z)dz + 0{e'^) 
for normalized initial data tpo G IIi/jqII = 1) via the quadrature formula 

[ {ao^t){z)n^„iz)dz = {1 +^) [ (ao$t)(z)H^o(2;)dz 

7lR2d J^2d 

d p 

j=l 

1 , d N N 

« —^ o $0(zfe) - — ^(a o $t)(wfc), 

k^l k^l 

where one samples the quadrature points according to the probability measures 

d 

Zi,..., zat ~ , Wi,..., Wat ~ i y] ; 

i=i 

see also for the sampling strategies used for the Hermite spectrograms. The 
rate of convergence for the above quadrature rule is proportional to for 

Monte Carlo samplings. For low discrepancy (Quasi-Monte Carlo) sampling the 
convergence is faster, that is, of the order log(7V)^‘^/iV. However, the literature 
on non-uniform Quasi-Monte Carlo sampling is scarce, and it seems to be an open 
question whether the transformed Halton sequences employed in our numerical 
experiments form indeed a low discrepancy set or just come very close to being so 
in practice, see also [T]. 

For the discretization of the Hamiltonian flow we apply the eighth-order 
symplectic splitting method from [271 Table 2.D], which is a suitable composition 
of the linear flows of 

.„d |« = " ,, . 

\p = 0 \p = -^V{q) 

Since our algorithm evolves an ensemble of classical trajectories, the use of sym¬ 
plectic time integrators is crucial; see also [H Fig. 4.2]. 


^All experiments have been performed with Matlab 8.3 on a 3.33 GHz Intel Xeon X5680 
processor. 
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Figure 3. Average errors (6.21 of the expectation values of vari¬ 
ous observables on the time interval [0, 20] for the new spectrogram 
method with initial Halton (left) and Monte Carlo (middle) sam¬ 
pling and results for the naive Husimi method with Halton sam¬ 
pling (right) for the torsional potential and Gaussian initial data 
centered at z = (1,0, 0,0). 


6.2. Two-dimensional torsional potential. Our first numerical experiments are 
conducted for the two-dimensional torsional potential 

V{qi,q 2 ) = 2 - cos(qi) - cos( 92 ), q € 

and different values of the semiclassical parameter e. As the initial state we consider 
the Gaussian wave packet -ipo = with phase space center z = (1,0, 0,0). This 
setup has already been considered in [giTiiiTiiTn]. We investigate the dynamics 
of the following symbols a : —>■ K, 

(1) Position: a{q,p) = qi and a{q,p) = q 2 , 

(2) Momentum: a{q,p) = Pi and a{q,p) = P 2 , 

(3) Kinetic and potential energy: a{q,p) = \\p\'^ and a{q,p) = V{q), 

(4) Total energy: a{q,p) = ijpp + V{q), 

and compare the outcome of the new algorithm with the naive, first-order Husimi 
approximation 


(6.1) (-0(t),op(a)V'(t)) = / {ao(^t){z)H^„(z)dz + 0{e). 

The left and middle panel of Figure confirm the second order accuracy of 
the new method for both Monte Carlo and Halton type samplings of the initial 
density p,^„. The right panel illustrates that the naive Husimi method is indeed 
only of first order in e. The time-averaged errors 


( 6 . 2 ) 



(V’ref(i),Op(a)V'ref(l)) 


1 -, d Af , N 

2 ^(ao$t)(zfe) + — ^(ao$t)(n;fc) 


N 




k^l 


dt 


are taken with respect to highly accurate grid based reference solutions tpref{t) 


^(t) of the Schrodinger equation; for details see Appendix B.l 
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The total energy error in Figurej^is smaller than the errors for other expectation 
values. Firstly, this can be explained by the fact that the total energy error is time 
independent, as explained in Remark |4.3[ and symplectic integrators are practically 
energy preserving on the time scales considered here. Secondly, the leading O(e^) 
term in the asymptotic expansion of the error 

{gz,op{h)g^)- h{w)gg^(w)dw = I - (1 - |A)'HgJ(w)dii; 

= [ h{w) (Wg^ - (1 - I A)(l + I A + ^A^)Wg) {w)dw + 0{e^) 

“ §2 / h{w)A'^yVg^{w)dw + 0{e^) 

JR^d 

can be bounded by the small constant 

1 ^ [ A^h{w)WgAw)dw\ < 

which follows from the special form of the torsional potential and the initial state. 
The total energy errors for small values of e are larger in the left panel of Figure 
than in the middle panel, since the relatively small number of Halton points gives 
rise to perceptible quadrature errors. 



Figure 4. The approximate dynamics of the potential ener¬ 
gies for the 32-dimensional Henon-Heiles system obtained by the 
Wigner, Spectrogram, and naive Husimi algorithms. 


6.3. Henon—Heiles dynamics for d — 32. Henon-Heiles type systems have 
been used for benchmark simulations with the multiconfiguration time-dependent 
Hartree method (MCTDH); see [52]. We follow the presentation in [Th] §5.B] by 
using the potential 

31 31 

V32(g) = + 1.8436 -b 

and the semiclassical parameter e = 0.0029, which is a model for the dynamics 
of a hydrogen atom on a high-dimensional potential energy surface that exhibits 
regions of chaotic motion. The quartic confinement guarantees that none of the 
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classical trajectories escapes to infinity. Moreover, as in [Sam], the initial state is 
a Gaussian wave packet = 9z centered at z = [q,p) with p = 0 and qj = 0.1215 
for all j = 1,..., 32. 

Since grid-based reference solutions are not available for this high-dimensional 
setting, we compare our method with the results obtained by the second order 
approximation 

= f {ao^t)iz)Wji,„iz)dz + 0{e‘^). 

Jm^d 


Note that, in this particular case, the initial Wigner function (z) is the Gauss¬ 
ian ( |5.1| ), and hence the Wigner algorithm does not pose a difficulty in the initial 
sampling. Figure shows the good agreement of the results from the Wigner and 
the spectrogram methods by means of the potential energy, and a considerable 
discrepancy with respect to the naive Husimi approximation (6.11. 


6.4. Escape from a cubic potential -well. Finally, we explore whether the new 
method is capable of describing the evolution of a quantum system that moves 
out of a potential well. For this purpose we consider the semiclassical Schrddinger 
operator H = — -|- E with the one-dimensional barrier potential 

V{q) = 2.328-q"^+q^+ 0.025q‘^, q G R, 

and e = 0.4642; see also Figure This Hamiltonian can be derived from the 
Schrddinger operator 

— -I- + O.la;^ 

with h = 1 from [33 IM] by applying the space rescaling x >->■ v'^OTa; and adding 
the confinement term 0.025 • q'^. The confinement prevents phase space trajectories 
from finite time blow up and guarantees that H is essentially self-adjoint. The 
global potential energy minimum H(xgiob) ~ —4765 is attained at a;giob ~ —28.4, 
and the confinement is very small in the region of interest close to the origin. 



Figure 5. The cubic barrier potential V with the barrier en¬ 
ergy 14 and the energy h{zo) of the trapped classical particle. 


As initial states we consider translated Hermite functions 

= fce{0,1,3,6} 

localized around zq = (0.4642,-1), which corresponds to the initial phase space 
center used in (33]. Since the associated classical energy h(zo) lies below the barrier 
energy 14 « 2.03, the classical particle is trapped in the well for all times. In 
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contrast, the quantum energy of the initial state Qzg is approximately 2.09, and the 
energies of the excited states are even higher. Consequently, the expected phase 
space center of the quantum particle will escape from the classical trapping region 
after short time. 

Figure [^displays the trajectories of the expected phase space centers obtained by 
the purely classical, the full quantum, and the semiclassical spectrogram dynamics 
for the four different initial states. The results from the spectrogram algorithm 
show decent qualitative agreement with the behavior of the quantum solution even 
though the semiclassical parameter e = 0.4642 is rather large. We also note that 
the results become more accurate for initial states of higher energy. 


Gaussian 
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Figure 6. Trajectories of the expected phase space centers ob¬ 
tained from quantum references and the new spectrogram algo¬ 
rithm for different initial states tpo = Tzg (pk ■ The dashed black line 
shows the periodic classical orbit associated with zq, and the solid 
black line illustrates the border of the trapping region. 


The spectrogram algorithm is also capable of describing the evolution of the 
probability that the quantum particle escapes the potential well and is found in the 
region (—oo, rcmax], where Xmax ~ —1.62 is the local maximum of the barrier poten¬ 
tial. To illustrate this property, we introduce the approximate escape probability 


P{t) = « IIV’tX(-oo,rr„„]lli2 
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with the smooth symbol r{q,p) = exp(-0.01/(g - a;max)^)X(-oo,a:„,,](<?), where XA 
denotes the characteristic function of the set A. P{t) can easily be approximated 
by the spectrogram algorithm, and accurate numerical references are available; see 
Appendix |B.3| Figure shows by means of two different initial states that the 

first excited state third excited state 




Figure 7. Approximate escape probabilities P{t) computed from 
a highly accurate numerical quantum reference, and results of the 
spectrogram and Wigner algorithms for the initial states Tz^ipi, 
and 


spectrogram algorithm yields a good qualitative picture of the evolution of escape 
probabilities. 

Appendix A. Sampling by the Gamma distribution 

A.l. Using the Gamma distribution. We consider the Hermite functions trans¬ 
lated by the Heisenberg-Weyl operator. 




z G 


i)2d 


Then, by the covariance property (2.6), the Hermite spectrogram takes the form 

)(w^) = we 

and by Lemma |5.2[ we only have to consider the two-dimensional probability den¬ 
sities 

2n 


hniWj - Zj) = 


27r£ • n! 


- Zj\^) 


with n G {/ci,..., fcd}. Specifically, we first translate by —Zj and then use the 
uniform distribution on [0, 27r] for the angular part combined with sampling from 

„2n+l 


r ::-: exp( — 


2"£"+in! 


for the radial part. Since 


fb ^2n+l 




exp(-i'^) dT 


2"e"+in!(2e)"+in! ' 

for all a, & > 0, we use the Gamma distribution with parameters n -|- 1 and 2e for 
the radial sampling. 
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In the particular case k = 0, (po = go and so is a Gaussian wave 

packet and 

Ms.M = (27re)"'^ (l + i - ik - exp(-^|u; - zf) , w G 

see §5.1| One can directly sample this 2d-dimensional probability density without 
factorizing its summands. After translation by —z, one decomposes into a radial 
part and an independent angular variable which is uniformly distributed on 
The radial density is given by 

"" (skd!®'’*"’'''"'’’ 

which corresponds to a Gamma distribution with parameters d + 1 and 2e. 

A.2. Monte Carlo sampling. Pseudorandom numbers uniformly distributed on 
a multi-dimensional unit sphere are obtained by sampling from multivariate normal 
distributions and subsequent normalization, while Gamma distributed pseudoran¬ 
dom samples only require the uniform distribution on the unit cube together with 
the (numerical) inverse of the cumulative Gamma distribution function. Hence a 
Monte-Carlo sampling of is straightforward. 

A.3. Quasi-Monte Carlo sampling. For the generation of Quasi-Monte Carlo 
points we have heuristically mimicked the Monte-Carlo procedure of §A.2| by re¬ 
placing the pseudorandom samples from the uniform distribution on the unit cube 
by a Halton sequence. Whether the resulting points are of low discrepancy with 
respect to the distribution pr^pk seems to be an open question not answered by 
the current literature, see e.g. [T]. 
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Appendix B. Numerical data 

B.l. Two-dimensional torsional system. Tablej^contains the number of initial 
sampling points and the computational time for the new spectrogram algorithm. 
In the case of Monte Carlo integration, the results in Figure are averaged over 
ten independent runs. For the time propagation we apply the above mentioned 
eighth-order symplectic integrator with time stepping 10“^. The parameters of the 
grid-based reference solver are collected in Table 


£ 

MC points 

comp, time 

Halton points 

comp, time 

10-1 

5-10^ 

23s 

5-10^ 

16s 

5•10-2 

3-103 

lm59s 

103 

33s 

10-2 

6-103 

7ml6s 

2 - 103 

lm59s 

5•10-3 

1.5 - 10® 

14ml5s 

8-103 

6m50s 

10-3 

10 - 10® 

68m30s 

2 - 10® 

18m31s 


Table 1 . Computational data for the execution of the new spec¬ 
trogram algorithm for the two dimensional torsional potential and 
Gaussian initial data on the time interval [0, 20]. The computation 
times are for one run only. They scale linearly with respect to the 
number of initial sampling points. 


6 

#timesteps 

comp, domain 

space grid 

10-1 

5 -103 

[-3,3] X [-3,3] 

1536 X 1536 

5 - 10-2 

5 - 103 

[-3,3] X [-3,3] 

1536 X 1536 

10-2 

7.5 -103 

[-2,2] X [-2,2] 

2048 X 2048 

5 -10-3 

10^ 

[-2,2] X [-2,2] 

2048 X 2048 

10-3 

10^ 

[-2,2] X [-2,2] 

2048 X 2048 


Table 2. Parameters of the grid-based reference solutions for the 
two-dimensional torsional potential. The discretization has been 
done by Fourier collocation in and Strang splitting in time. 


B.2. Henon—Heiles system for d — 32. For the initial sampling we used 2^^ 
Halton type points for all three initial densities, that is, the Wigner and the Husimi 
functions and the new density /i^g. The time stepping of the eighth order time 
integrator is 2 • 10“^. 

B.3. One-dimensional cubic well. For the spectrogram algorithm we employed 
2^^ Halton points and the eighth-order integrator with time stepping 10“^, which 
results in a computational time of 2 seconds. The quantum references are generated 
by means of a Strang splitting with 2^^ Fourier modes on the interval [—40,4]. 
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